Crank-Nicolson in 2D

See Sections 9.8 (LOD) & 9.8.2 (ADI) for other methods.

In this notebook we will develop numerical methods to solve

\begin{align*} \begin{cases} u_t(x,y,t) = u_{xx}(x,y,t) + u_{yy}(x,y,t)\\ u(x,y,0) = \eta(x,y),\\ u(0,y,t) = g_0(y,t),\\ u(1,y,t) = g_1(y,t),\\ u(x,0,t) = h_0(x,t),\\ u(x,1,t) = h_1(x,t). \end{cases} \end{align*}

Using the method of lines, we obtain for $1 \leq i,j \leq m$, $h = 1/(m+1)$,

\begin{align*} U'_{i,j}(t) = \frac{ U_{i-1,j}(t) + U_{i+1,j}(t) + U_{i,j-1}(t) + U_{i,j+1}(t) - 4 U_{i,j}(t)}{h^2}. \end{align*}

Applying the trapezoid method to this we obtain 2D Crank-Nicolson. If we use forward Euler we obtain (clearly) an explicit method. The von Neumann stability analysis indicated that this explicit method is stable if

\begin{align*} \frac{k}{h^2} \leq \frac{1}{4}. \end{align*}

We will confirm this to be the case. If the grid points $(x_i,y_j)$ are arranged by the ordering $\prec$ defined by

\begin{align*} (x_{i'},y_{j'}) \prec (x_{i},y_{j}) \quad\text{if}\quad j' < j, \quad \text{and},\\ (x_{i'},y_{j}) \prec (x_{i},y_{j}) \quad\text{if}\quad i' < i. \end{align*}

This orders first by the $y$ values and then for a given $y$ value it orders by increasing $x$ value. This is just the ordering given in Figure 3.2(a). With this ordering the discrete operator from the MOL approach is given by (3.12):

\begin{align*} U'(t) = \frac{1}{h^2}A U(t),\quad \text{where} \quad A= \begin{bmatrix} T & I\\ I & T & I \\ & I & T & I\\ && \ddots & \ddots & \ddots \\ &&& I & T \end{bmatrix}, \quad T = \begin{bmatrix} -4 & 1\\ 1 & -4 & 1 \\ & 1 & -4 & 1\\ && \ddots & \ddots & \ddots \\ &&& 1 & -4 \end{bmatrix}. \end{align*}

Note that $A$ is an $m^2 \times m^2$ matrix. In our calculations we will have $r = k/h^2$ for the explicit method and $r = k/(2h^2)$ for the implicit method. With that in mind, we have two important matrices

\begin{align*} A_l = I - r A, \quad A_r = I + r A. \end{align*}

And $A_l$ is only used in the implicit methods.

The building blocks

  1. Build the grid
  1. Convert values on the grid to vectors
  1. Handle the boundary conditions

For example, grid points $1,2,\ldots,m$ are connected to $h_0(x,t)$ and grid points $1,m+1,2m+1,\ldots$ are connected to $g_0(y,t)$.

  1. Construct $A$
  1. Initial/boundary data and plotting helper

An actual run

Test stablity with explicit method

The implicit method

We need an implementation of the conjugate gradient algorithm. We could be using the locally one dimensional (LOD) or the alternating direction implicit (ADI) methods. But I wanted to show you how one would incorporate an iterative method for sparse systems into a PDE solver.

Preconditioning?

Timings

Create a gif with a high framerate